
;;================================================================
;; Routine to generate the ATLASGAL - BGPS latitude comparison for the
;;   overlap between the two surveys.  BOTH catalogs have been
;;   restricted to the region -10 < l < 21 (something not done in
;;   Fig. 9 of {Contreras:2013}).
;;
;; This routine makes flux-weighted distributions to remove the
;;   vagaries of catalog algorithm junk.
;;


;; BGPS & ATLASGAL data...
csc = read_mrt('./ancillary/ATLASGAL_csc.dat')
s = read_bgps_csv(/ver)
s.flux *= 1.5

;; Cut each catalog to the overlap region -10.5 >= l >= 21
ind = where(csc.glon GE 100)
csc[ind].glon -= 360.
ind = where(csc.glon GE -10.5)
csc = csc[ind]

ind = where(s.glon_peak GE 250)
s[ind].glon_peak -= 360.
ind = where(s.glon_peak LE 21)
s = s[ind]

print,mean(s.glat_peak),total(s.glat_peak*s.flux)/total(s.flux)
print,mean(csc.glat),total(csc.glat*csc.fint)/total(csc.fint)


nbin = 30
xp = dindgen(nbin)/10.d - 1.45d
print,m4_stat(xp)

print,n_elements(csc),n_elements(s)

fa = dblarr(nbin)
fb = dblarr(nbin)


FOR ii=0,nbin-1 DO BEGIN
   
   ia = where( abs(csc.glat - xp[ii]) LE 0.05, na)
   ib = where( abs(s.glat_peak - xp[ii]) LE 0.05, nb)
   print,xp[ii],na,nb
   
   fa[ii] = na GT 0? total(csc[ia].fint) : 0.d
   fb[ii] = nb GT 0 ? total(s[ib].flux) : 0.d
   
ENDFOR

myps,'./emaf_paper/plots/atlasgal_fluxglat.eps'
deg = cgSymbol('deg')

fa /= 1.d3
fb /= 1.d3
fb *= 2.25454; * 3.3851

cgPlot,charsize=1.0,xp,fa,psym=10,xr=[-1.55,1.55],/xst,thick=5,/nodata,$
       xtit='Galactic Latitude [deg]',ytit='Integrated Source Flux Density per 0.1 deg bin [kJy]'

FOR ii=0,nbin-1 DO $
   cgPolygon,xp[ii]+[-0.05,-0.05,0.05,0.05],[0,fa[ii],fa[ii],0],$
             color='BLK4',/fill,fcolor='BLK4',thick=5

scale = 1.d; * max(fa) / max(fb)


yfa = mpfitpeak(xp,fa,Aa,nt=3)
yfb = mpfitpeak(xp,fb,Ab,nt=3)
;; scale = Aa[0]/Ab[0]
print,'SCALE: ',scale

xc = findgen(4001)/1000. - 2.
yca = gauss_1(xc,Aa)
ycb = gauss_1(xc,Ab)

cgOplot,xp,fb*scale,psym=10,color='Black',thick=5
cgOplot,xc,yca,thick=5,linestyle=2,color='blk5'
cgOplot,xc,ycb*scale,thick=5,linestyle=0,color='black'


cgAxis,xaxis=0,/xst,xtickformat='blank_axis'

;; ATLASGAL Legend
al_legend,/top,/left,box=0,textcolor=['Opposite','BLK4','BLK6'],$
          ['ATLASGAL','CENTROID = '+string(Aa[1],format="(F0.2)")+deg,$
           'FWHM = '+string(Aa[2]*2.355,format="(F0.2)")+deg],linsize=0.5,$
          color=['Background','BLK4','BLK6'],thick=5,linestyle=[0,0,2]

;; BGPS Legend
al_legend,/top,/right,box=0,$
          textcolor=['Opposite','Opposite','Opposite'],$
          ['BGPS','CENTROID = '+string(Ab[1],format="(F0.2)")+deg,$
           'FWHM = '+string(Ab[2]*2.355,format="(F0.2)")+deg],linsize=0.5,$
          color=['Background','Opposite','Opposite'],thick=5,linestyle=[0,0,0]

myps,/done

END
